function smp = srf_CurvatureMap(hfile)
% SRF::CurvatureMap  - compute a curvature map for a surface
%
% FORMAT:       smp = srf.CurvatureMap;
%
% No input fields.
%
% Output fields:
%
%       smp         SMP with one map containing the curvature info

% Version:  v0.7f
% Build:    8110521
% Date:     Nov-05 2008, 9:00 PM CET
% Author:   Jochen Weber, SCAN Unit, Columbia University, NYC, NY, USA
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% check arguments
if nargin ~= 1 || ...
    numel(hfile) ~= 1 || ...
   ~isBVQXfile(hfile, 'srf')
    error( ...
        'BVQXfile:BadArguments', ...
        'Invalid call to %s.', ...
        mfilename ...
    );
end
sc = bvqxfile_getscont(hfile.L);
bc = sc.C;

% prepare arrays for computation
n = bc.NrOfVertices;
t = bc.TriangleVertex;
vo = double(bc.VertexCoordinate(:, :));
vn = vo - 0.01 * double(bc.VertexNormal(:, :));

% compute area per triangle
tvo = cat(3, vo(t(:, 1), :), vo(t(:, 2), :), vo(t(:, 3), :));
tvo = cat(3, ...
    tvo(:, :, 2) - tvo(:, :, 1), ...
    tvo(:, :, 3) - tvo(:, :, 2), ...
    tvo(:, :, 1) - tvo(:, :, 3));
tvo = squeeze(sqrt(sum(tvo .* tvo, 2)));
so = sum(tvo, 2) ./ 2;
ao = sqrt(so .* (so - tvo(:, 1)) .* (so - tvo(:, 2)) .* (so - tvo(:, 3)));

% the same at 0.01 along normal
tvn = cat(3, vn(t(:, 1), :), vn(t(:, 2), :), vn(t(:, 3), :));
tvn = cat(3, ...
    tvn(:, :, 2) - tvn(:, :, 1), ...
    tvn(:, :, 3) - tvn(:, :, 2), ...
    tvn(:, :, 1) - tvn(:, :, 3));
tvn = squeeze(sqrt(sum(tvn .* tvn, 2)));
sn = sum(tvn, 2) ./ 2;
an = sqrt(sn .* (sn - tvn(:, 1)) .* (sn - tvn(:, 2)) .* (sn - tvn(:, 3)));

% build curvature index (per triangle)
crv = 1 - an ./ ao;

% average for each point of triangle
vv = zeros(n, 1);
vd = zeros(n, 1);
for vc = 1:3
    for tc = 1:size(t, 1)
        tr = t(tc, vc);
        vv(tr) = vv(tr) + crv(tc);
        vd(tr) = vd(tr) + 1;
    end
end

% average correctly
vv = vv ./ vd;

% create smp
smp = bless(BVQXfile('new:smp'), 1);
smc = bvqxfile_getcont(smp.L);
smc.NrOfVertices = n;
smc.NameOfOriginalSRF = sc.F;
smc.Map(1).DF1 = 0;
smc.Map(1).RGBLowerThreshPos = [0, 85, 170];
smc.Map(1).RGBUpperThreshPos = [64, 128, 255];
smc.Map(1).RGBLowerThreshNeg = [0, 170, 85];
smc.Map(1).RGBUpperThreshNeg = [64, 255, 128];
smc.Map(1).LowerThreshold = 0;
smc.Map(1).UpperThreshold = 1.5;
smc.Map(1).Name = 'Curvature';
smc.Map(1).SMPData = 100 * vv;
bvqxfile_setcont(smp.L, smc);
